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An old and apparently persistent problem in numerical lattice QCD is that the simulations tend to 
get trapped in a sector of fixed topological charge when the lattice spacing is taken to zero. The 
effect sets in very rapidly and may invalidate the simulation results in certain cases. In this talk, 
the issue is discussed using the Wilson flow as a tool. The flow has a simple scaling behaviour 
and allows one to understand how exactly the topological sectors emerge in the continuum limit. 
Further studies however suggest that the observed slowdown of the simulations at small lattice 
spacings is only partly caused by the emergence of the sectors. 
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1. Introduction 

At present all widely used simulation algorithms for lattice QCD with two or more flavours of 
sea quarks are based on some version of the HMC algorithm []l|] (see ref. [^] for a recent review). 
The computer time required for such simulations depends on the time needed for an update step, 
the integrated autocorrelation times of the observables of interest and the desired statistics. As a 
function of the lattice spacing a, and if the algorithm is implemented as usual, the first factor scales 
roughly like cT 5 , but the autocorrelation times tend to be difficult to determine and are often not 
known reliably if they are not small. 

Large autocorrelation times are typically observed in the case of quantities related to the topo- 
logical charge of the gauge field. While the definition of the charge on the lattice is ambiguous to 
some extent, experience suggests that the autocorrelation times are largely insensitive to the exact 
choices one makes. For illustration a time series of charge measurements is plotted in fig. [T[ The 
presence of strong autocorrelations is evident in this example, particularly so in the shaded region, 
where the charge never changes sign and instead oscillates around a value of —10. Note that most 
QCD simulations published to date are much shorter than the run shown in the figure. 

The integrated autocorrelation times of physical quantities tend to grow when the continuum 
limit is approached. Depending on the theory, the algorithm and the observable, the asymptotic 
scaling behaviour can be very different, but is normally power-like in the inverse lattice spacing 
with an exponent z < 2. Del Debbio, Panagopoulos and Vicari [||] however found some time ago 
that the topological charge in the SU(3) gauge theory has a much more rapidly (perhaps expo- 
nentially) increasing autocorrelation time if the standard link-update algorithms are used. Further 
studies by Schaefer, Sommer and Virotta [Q] later showed that the HMC algorithm is similarly in- 
efficient and there is ample evidence that the problem persists when the sea quarks are added to the 
theory @-@|. 

In this talk, a few steps are taken towards a better understanding of the dynamics of the HMC 
algorithm and the mechanism that leads to the dramatic slowdown of the simulations at small lattice 
spacings. Issues to be addressed are how exactly the topological sectors emerge in the continuum 
limit, what the slow modes of the gauge field might be and whether perhaps there is a simple way 
out. The discussion is largely based on a new tool, the Wilson flow [[7], ^J, which is of some interest 
in its own right. 

2. Autocorrelations 

2.1 Scaling behaviour 

Systematic scaling studies of autocorrelation times in lattice QCD have so far been limited 
to the case of the pure gauge theory [^, ||]. In all these studies, the autocorrelation time of the 
topological charge turned out to increase at least like a~ 5 at lattice spacings a below 0. 1 fm or so. 
The algorithms considered (the HMC, the DD-HMC and the well-known link-update algorithms) 
appear to behave similarly in this respect. In particular, once the link- visiting frequency is divided 
out, the block size used in the DD-HMC algorithm does not have a significant influence on the 
autocorrelation times [Q]. 
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Figure 1: History of the topological charge in a simulation of the pure SU(3) gauge theory on a 64 x 32 3 
lattice with spacing a = 0.07 fm. The plain HMC algorithm was used in this test run, with a trajectory length 
of 2 units of molecular-dynamics time and an acceptance rate of 83%. 

The onset of the rapid growth of the autocorrelation time of the topological charge, and thus 
its value at a given lattice spacing, however depends on the chosen lattice action and the simulation 
algorithm. Increasing the length of the molecular-dynamics trajectories in the HMC algorithm can 
be beneficial for example, while the addition of six-link terms to the Wilson plaquette action 
may have an adverse effect [Q, ^]. 

When the sea quarks are included in the simulations, the situation becomes considerably more 
complicated, because the autocorrelations may now also depend on the number of quark flavours, 
the quark masses and chosen the fermion action These dependencies and the one on the 

lattice spacing remain to be studied in detail, but the experience made so far (which is sometimes 
only based on a visual inspection of measurement histories) shows that the autocorrelation time of 
the topological charge is again very rapidly growing when the lattice spacing is reduced from 0. 1 
fm to 0.05 fm or even smaller values. 

If the leading exponential autocorrelation time is assumed to grow proportionally to a~ 5 , the 
total computational effort required for HMC simulations of QCD is expected to scale like a -10 at 
fixed physics, i.e. the cost of the simulations increases by about three orders of magnitude when 
the lattice spacing is divided by 2. This estimate may be a bit pessimistic, but it is quite clear that 
the simulations required for safe extrapolations to the continuum limit are extremely challenging. 
Algorithmic improvements or viable ways of bypassing the slowing down of the simulations are 
certainly highly desirable at this point. 

2.2 Autocorrelation effects in short runs 

A question often asked in this context is whether the simulations really need to be very much 
longer than the leading exponential autocorrelation time. In particular, if the quantities of interest 
are only weakly coupled to the slow modes of the algorithm, the results obtained in shorter runs 
may conceivably be correct within statistical errors. 

In general, the expectation values calculated in such short runs must be expected to be biased 
to some extent. Considering again the measurement history plotted in fig. [T], for example, it is clear 
that the topological charge is incorrectly sampled in this case if runs not very much longer than 5000 
molecular-dynamics time units are performed. The expectation values of most observables are then 
affected by terms inversely proportional to the space-time volume []9[ [Kj. Physical quantities like 
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Figure 2: Example of a normalized autocorrelation function p(f) that is dominated by two modes with 
exponential autocorrelation times To 3> T\ . 

the masses of the rj and rj' mesons, but also hadronic matrix elements of pseudo-scalar densities 
are likely to be strongly sensitive to these effects. 

When arguing for short runs, one would need to provide a practical procedure that allows the 
effects of the slow modes on the calculated expectation values to be estimated. Whether they are 
in fact negligible in the cases of interest is otherwise difficult to tell and one is left with results that 
may or may not be correct. 

2.3 Statistical error estimation in presence of slow modes 

Large exponential autocorrelation times can also lead to an underestimation of integrated au- 
tocorrelation times and thus of the associated statistical errors. To illustrate this point, consider a 
normalized autocorrelation function 

p(t) = |c | 2 e- f / To + |d 1 2 e-'/ T ' + . . . (2.1) 

which is dominated by two eigenmodes of the simulation transition probability, a slow mode with 
exponential autocorrelation time To and a fast mode with autocorrelation time Ti <C To (see fig. |2|). 
The coefficients |co| 2 and |ci| 2 in this formula measure how strongly these modes couple to the 
observable considered. Note that the integrated autocorrelation time 

/ Tint-To|c | 2 + Ti|ci| 2 + ... (2.2) 

is equal to the area under the curve shown in the figure. Whether Tj nt can be easily estimated now 
depends on whether the area under the tail of the curve dominates or not. If it does not, i.e. if 
^okol 2 ^ T il c i| 2 ' the relevant exponential autocorrelation time is Ti and the statistical error is 
correctly obtained from simulations a few hundred times longer than Ti . Runs very much longer 
than To are however required in the other case to be able to control the situation. 

Note that the contribution of the slow mode to the statistical error is not guaranteed to be 
negligible even if the coupling |co| 2 is very small. These cases are actually particularly difficult to 
treat correctly, because the autocorrelation function in the tail is very small, while the area under 
the tail may not be so. Since the autocorrelation function itself can only be calculated up to some 
statistical uncertainty, it may then be quite impossible to exclude this case, except when To is known 
or can at least be bounded from above [|J]. 
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Figure 3: Equation (3.1) generates a flow in the space of lattice gauge fields. On a finite lattice, the existence, 
smoothness and uniqueness of the flow is rigorously guaranteed at all flow times Moreover, at any given 
time t, the mapping U Vt is an invertible transformation of field space, because the flow can be integrated 
both forward and backward in time. 



3. Wilson flow 



The high-frequency components of the gauge field are usually weakly coupled to the slow 
modes of the HMC algorithm and tend to be efficiently updated. What exactly one means by the 
smooth components of the field is not entirely clear, however, and even less so which of their 
properties (apart from the topological charge) are slowly sampled by the algorithm. As explained 
in the following, the Wilson flow ^ ^] is a renormalizable smoothing operation that allows these 
questions to be addressed on theoretically solid ground. 



3.1 Flow equation 

For any given lattice gauge field U(x,n), the first-order differential equation 

V t ( X ,n) = -gl{d X!fl S w (V t )}V t ( X ,n), V t (x,n)\ t=0 = U( X ,Ll), (3.1) 

defines a trajectory V t (x,n) of fields parameterized by the "flow time" t (see fig. [|; differentiation 
with respect to t is abbreviated by a dot). In this equation, S w (V t ) denotes the Wilson plaquette 



action Jl lfl of the field V t at gauge coupling go and d xll S w (V t ) its (Lie algebra valued) variation 
with respect to the link variable V t (x,jU). Note that the coupling cancels in the flow equation. 

Along the Wilson flow, the plaquette action decreases monotonically, 5 W < 0, and the gauge 
field tends to become smoother. The flow is in fact generated by infinitesimal "stout" link-smearing 



steps [ ]12| ] and thus shares some properties with this popular smoothing procedure. Eventually 
the flow drives the field towards the stationary points of the action, but contrary to what may be 
assumed, the large-time regime is a highly non-perturbative one in QCD. In particular, it may not 
be meaningful to study the large-time behaviour of the flow separately from the continuum limit of 
the theory. 
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3.2 QED — a soluble case 

In the continuum limit, the flow equation ( |3~l| ) assumes the form 

B jl =D v G v ^ l , B fl \ t=Q =A fl , (3.2) 

where A^(x) is the fundamental gauge field, B^(t,x) the time-dependent gauge field and G^ v 
the associated field tensor. In QED the equation is a linear diffusion equation, whose solution is 
given by 

B^(t,x) = J d 4 yK t {x-y)A fl (y)+gmgt terms, K,(z) = (3.3) 

The flow thus averages the gauge field over a spherical range with mean-square radius equal to y/St. 
Note that, as is already clear from eq. ( |3.2| ), the flow time t has engineering dimension [length] 2 . 

Since the field generated by the flow is linearly related to the fundamental gauge field, its 
correlation functions 

(B lll (t,xi)...B IJn (t,x„)) = e oJ A A yi...d 4 y n K 1 (x l -yi)...K t (x n -y n ) 

xG (ji,. . . ,y, 1 ) A(l ... A( , I + gauge terms (3.4) 

are proportional to the (full, bare) photon rc-point functions Go(yi,. . . ,3 ; ,i)/j 1 ...^„- A power of the 
bare electron charge eo appears in this equation, because the canonically normalized photon field 
is eoAn rather than (which is normalized so that the covariant derivatives do not involve the 
charge). Note also that the heat kernels K t (xk —yk) play the role of smooth test functions in this 
formula. In particular, the renormalization of the correlation function is achieved simply by renor- 
malizing the bare charge and the photon «-point function according to 

e = Z 3 i/2 e R , G = Z'^ 2 G R . (3.5) 

The fact that both eo and Go renormalize with the same renormalization constant Z3 is a conse- 
quence of the gauge Ward identity in this theory (a gauge-invariant regularization is assumed here). 
Since e^Go = e R G R , this shows that the field Bn(t,x) is, at all positive flow times t and up to its 
gauge degrees of freedom, a renonnalized smooth gauge field. 

Provided the bare charge is expressed through the renormalized one, the correlation functions 
of the field tensor G^ v thus do not require any renormalization and converge to well-defined smooth 
functions of the space-time coordinates when the regularization of the theory is removed. More- 
over, since QED is asymptotically free at low energies, the behaviour of the correlation functions 
at large flow times is described by leading-order perturbation theory. A short calculation then leads 
to the formula 

lim{? 2 (G AIV G Alv )} = ^, (3.6) 

which shows that the field obtained by the Wilson flow contains some interesting physical infor- 
mation. 
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Figure 4: Plot of t 2 (E) in the pure SU(3) gauge theory as a function of the flow time t given in units of 
the Sommer scale [|l3j p4[ ]. The thick solid line was obtained through numerical simulation of a 96 x 48 3 
lattice with spacing a = 0.05 fm. One-loop perturbation theory [eq. (3.8)] and the known value of the A- 
parameter jl5| ] yield the dashed curve with an error margin (grey area) deriving from the error on rgA. 



3.3 Properties of the Wilson flow in QCD 

In QCD the flow equation is non-linear and the renormalization of the theory is much more 
complicated than in QED. Whether the Wilson flow generates a renormalized gauge field is how- 
ever a question that can be studied in perturbation theory. In particular, the expectation value of the 
gauge-invariant density 

E = -itr{G MV G MV } (3.7) 

can be easily worked out to next-to-leading order in the gauge coupling. Dimensional regularization 
may be used in this calculation and it then turns out that (E) does not require renormalization to 
this order, i.e. the divergent terms are all canceled by the renormalization of the coupling. In the 
MS scheme, the one-loop formula obtained in this way is [JSj] 

(E) = Jt- a ( q ){l+k ia (q) + ...}, <7 = (8*r 1/2 , (3-8) 
Jti = 1 .0978 + 0.0075 x Nf, (3.9) 

where Nf denotes the number of massless sea quarks and the running coupling a(q) is evaluated at 
a momentum scale q equal to the inverse of the leading-order smoothing range of the flow. 

Beyond perturbation theory, the expectation value of E can be computed straightforwardly 
using the lattice formulation of the theory and numerical simulations. An accurate numerical inte- 



gration of the flow equation (3.1) is required in these calculations, but the computer time needed for 
the integration is negligible in practice if a suitable higher-order integrator is used The result 
of a computation along these lines is plotted in fig. [| in a range of the flow time corresponding to 
smoothing ranges from about 0.2 to 0.5 fm (the statistical errors are not visible on the scale of the 
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Figure 5: Scaling behaviour of the reference time to in the pure SU(3) gauge theory. The data (full points) 
were obtained on lattices of size 48 x 24 3 , 64 x 32 3 and 96 x 48 3 with spacings approximately equal to 0.1, 
0.07 and 0.05 fm, respectively. The continuum limit (open point) is reached through a linear fit in a 2 of the 
three data points (solid line). 



plot). Since QCD is asymptotically free, the perturbation expansion ( [3lj| ) only applies at small flow 
times, but as can be seen from the curves shown in fig. ||, the transition from the small-time to the 
non-perturbative regime is very smooth. Note that t 2 {E) increases roughly linearly with t in the 
non-perturbative regime, at least so in the range shown in the plot, a behaviour which is completely 
different from the one in QED. 

If (E) does not require renormalization, one expects t 2 {E) to be a universal function of tjr\ up 
to lattice effects that vanish proportionally to a positive power of the lattice spacing. The scaling 
can be checked by introducing a reference time to through the implicit equation 



0.3 



(3.10) 



(see fig. ^]). Simulations at three values of the lattice spacings then show that the dimensionless 
ratio v^S/o/ro does in fact smoothly converge to the continuum limit, the lattice-spacing effects 
being less than a percent in the range covered by the simulations (fig. |5|). Little doubt thus remains 
that the Wilson flow maps the gauge field to a renormalized smooth field as in QED. 



4. Topological sectors 

The question of how exactly the topological sectors emerge in lattice QCD can now be an- 
swered by performing the field transformation U — > V = V tQ in the functional integral, to being the 
reference flow time introduced in the previous section. The smoothing range is about 0.5 fm in 
this case, i.e. on average the fluctuations of the gauge field with wavelengths up to the confinement 
radius are smoothed out. Somewhat surprisingly, the Jacobian of the transformation can be analyti- 
cally expressed through the Wilson action [JJ. The expectation value of any observable ff(U) then 
assumes the form 

(@) = ^jw\m)^ iy \ (4.i) 
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Figure 6: Probability for the value of s p [eq. (4.3)] on a given plaquette p to be above a specified threshold s. 
The three curves were obtained in the pure SU(3) gauge theory with Wilson plaquette action, using the same 
representative ensembles of gauge-field configurations as in the scaling test reported in fig. 5. 



where the action is given by 



§(V) = S(U)+ l 44 r°d^ w (V f ). (4.2) 
3a z Jo 



In this formula, the fields U and V t are considered to be functions of V and S(U) denotes the action 
of the theory (including the quark determinants if any) before the transformation. 

Both terms in the action ( [Q| ) tend to suppress large values of the plaquette observable 

s p = Reti{\-V(p)} (4.3) 

(where V(p) is the product of the link variables V(x,fi) around the plaquette p) and thus force the 
plaquette loops to be close to unity. Since s p is a field of dimension 4, representing the square of 
the gauge-field tensor in the plaquette plane, its expectation value is in fact expected to scale like 
a 4 in the continuum limit. Numerical studies confirm this behaviour and show that large values of 
s p are indeed very strongly suppressed (see fig. ||). 



Many years ago, a theorem was established [16, 17] which states that the space of all lattice 
gauge fields satisfying a certain smoothness condition decomposes into topological charge sectors 
very much like the space of continuous fields in the classical continuum theory. In particular, the 
theorem applies to the subspace of SU(3) gauge fields V which satisfy 

s p < 0.067 for all plaquettes p. (4.4) 

According to the empirical results reported in fig. ^, and since the number of plaquettes in a fixed 
physical volume grows proportionally to aT A , the fields that do not fulfill this condition occur with 
a probability of order a 6 . In the functional integral (fO|), the weight of these configurations (and 
thus of the region of field space, where the assignment of the topological charge is ambiguous) is 
therefore rapidly decreasing when the lattice spacing is taken to zero. In lattice gauge theory, the 
emergence of the topological sectors is thus seen to be a dynamical phenomenon, which sets in 
close to the continuum limit. 
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Table 1: Integrated autocorrelation times in the pure SU(3) gauge theory with Wilson action, measured on a 
48 x 24 3 lattice with spacing a = 0.1 fm (left table) and on a 64 x 32 3 lattice with spacing a = 0.07 fm (right 
table). In both cases, the HMC algorithm with a trajectory length of 2 and an acceptance rate of 83% was 
used. All autocorrelation times and the trajectory length are given in units of molecular-dynamics time. 

5. Slow modes of the HMC algorithm 

The field transformation considered in the previous section maps the molecular-dynamics tra- 
jectories generated by the HMC algorithm to trajectories of the transformed field V = V t(t . At lattice 
spacings where the regions of field space "between the topological sectors" are strongly suppressed, 
it is unlikely that such a trajectory leads from one sector to another, because the fields along the 
trajectories are, to a good approximation, distributed according to their weight in the functional in- 
tegral. The emergence of the sectors and the freezing of the topological charge in HMC simulations 
are thus directly related to each other. 

5.1 Do the smoothed fields move slowly? 

The algorithm may however slow down for other reasons as well and there is no guarantee that 
the gauge field is efficiently updated in a fixed sector. In particular, the smooth fields obtained by 
the Wilson flow at flow time t > may conceivably be slowly moved through configuration space, 
in which case the autocorrelation times of local observables constructed from these fields (such as 
E) are expected to be large. 

The measured autocorrelation times of E listed in table 1 are actually monotonically increasing 
with the flow time t and reach values up to 10 times the autocorrelation time of the usual plaquette 
observable at t = 0. On the other hand, when the lattice spacing is reduced, they do not appear to 
grow as rapidly as the autocorrelation times of the topological charge Q (in these studies, E and 
the charge density were both defined using the same symmetric plaquette-loop expression for the 
gauge-field tensor at the specified flow time). The relatively slow updating of the smoothed field 
and the charge freezing are thus two different effects, where one is not obviously driving the other. 

5.2 Open boundary conditions 

One may be inclined to conclude at this point that the emergence of the topological sectors at 
small lattice spacings is the principal and perhaps only cause for the presence of very rapidly grow- 
ing autocorrelation times in HMC simulations. However, when choosing open instead of periodic 
boundary conditions in the physical time direction, the topological sectors disappear and the space 
of smooth fields becomes connected. One thus expects to observe only moderately increasing au- 
tocorrelation times in this case if the slowdown of the algorithm is indeed mainly caused by the 
separation of the sectors. 
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periodic 68(6) 
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periodic 615(90) 
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Table 2: Comparison of autocorrelation times measured on lattices with periodic and open boundary condi- 
tions in the physical time direction. The lattices are the same as in table 1 (left: a = 0.1 fm, right: a = 0.07 
fm), the HMC parameters are also the same and the flow time t is equal to 0.7 x to in all cases. 

A lattice with time extent T and open boundary conditions in time does not wrap around in 
this direction, i.e. there are no terms in the action which couple the field variables at time xq = to 
those at the largest time xq = T — a. In a simulation program for periodic lattices, these boundary 
conditions can often be implemented simply by setting the time-like link variables at the largest 
time to zero. Note that open boundary conditions preserve the gauge symmetry. In the continuum 
theory, they amount to imposing Neuman boundary conditions, 

Fok(x)\ Xo=0 = F ok {x)\ XQ=T = 0, k = 1,2,3, (5.1) 

on the gauge field and homogeneous Dirichlet (Schrodinger functional) boundary conditions on the 
quark fields, F^ v being the field tensor of the fundamental gauge field. The field space becomes 
topologically trivial when these boundary conditions are chosen. In particular, instantons can be 
smoothly moved in and out of the volume through the boundaries at xo = 0, T. 

The figures listed in table || show that the situation is in fact somewhat improved when passing 
from periodic to open boundary conditions (in both cases, Q denotes the sum of the topological 
charge density over all lattice points). As a function of the lattice spacing, the autocorrelation time 
of Q is however rapidly growing independently of the chosen boundary condition. On the periodic 
lattices, the quantization of the topological charge thus appears to be only one of several causes of 
the slowdown of the HMC algorithm. 

5.3 Instantons and the chiral limit 

One of the distinguishing features of non-abelian gauge theories is the existence of the instan- 
ton solutions of the classical field equations. The fact that the action has a large manifold of nearly 
stationary points (the configurations built from many distant instantons and anti-instantons) may 
conceivably play an important dynamical role in these theories, but so far it proved to be difficult 
to actually show this. In the present context, the observation is perhaps of some relevance, because 
the fields generated by the Wilson flow at large flow times tend to be close to the stationary points 
of the action. 

In QCD with light sea quarks, the fluctuations of the topological charge are suppressed and the 
smooth fields obtained by the Wilson flow may therefore be characteristically different from the 
ones in the pure gauge theory. The slow modes of the HMC algorithm consequently need not be the 
same and the autocorrelation times of the usual quantities may turn out to have a significant depen- 
dence on the quark masses. However, as already mentioned, systematic studies of autocorrelation 
times await to be performed in QCD with dynamical quarks. 
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6. Conclusions 

Numerical lattice QCD rests on the assumption that simulations at sufficiently small lattice 
spacings can eventually be performed to be able to control the discretization errors. It is therefore 
of central importance to overcome the poor scaling behaviour of the currently available simulation 
algorithms. Finding better algorithms is difficult, however, and probably requires a more detailed 
understanding of why the HMC algorithm slows down when the lattice spacing is taken to zero. 

The Wilson flow is a useful tool in this context, because it allows the properties of the gauge 
field at different length scales to be studied consistently with the renormalization of the theory. In 
particular, the way in which the topological sectors emerge close to the continuum limit is made 
transparent through the flow. The emergence of these sectors is certainly one of the principal causes 
of the slowdown of the HMC algorithm, but the fact that the situation improves only slightly on 
lattices with open boundary conditions is a bit surprising and remains unexplained. It seems safe 
to conclude, however, that there are important further sources of inefficiency and that the problem 
is likely to persist at fixed topological charge (the rapid slowdown of the algorithm is in this case 
expected to be revealed when considering the sum of the charge density over half the lattice). 

Choosing open boundary conditions in the physical time direction is nevertheless an inter- 
esting option, because the barriers between the sectors disappear in this case, while the transfer 
matrix (and thus much of the physics) is the same as with periodic boundary conditions. From 
the point of view of the simulation algorithms, the absence of the barriers represents an important 
simplification. Tunneling transitions are then not required anymore and one is left with the task of 
finding an algorithm that moves the field V t (at, say, t = to) efficiently through the space of smooth 
configurations. 
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